#include "bc_dynamicValue.hpp"
#include "meshTypes.hpp"
#include "field.hpp"
#include "constantField.hpp"
#include "fieldFunctions.hpp"
#include "boundaryPatch.hpp"
#include "simulationTime.hpp"
#include <iostream>

namespace hms
{


void assertSortedValues(
	[[maybe_unused]] const TimeSeries& values
){
	assert( std::is_sorted(
		values.cbegin(),
		values.cend(),
		[]( const auto& pair1, const auto& pair2 ){
			return pair1.first < pair2.first;
		}
	) );
}

#define CHEAD(RETURN_TYPE) \
template<typename MeshType> \
RETURN_TYPE DynamicValue<MeshType>::

CHEAD() DynamicValue(
	Field<MeshType>& field,
	const BoundaryPatch& patch,
	const Time& timeArg
) :
	BaseType( Type::dynamicValue, field, patch )
{
	time(timeArg);
}

CHEAD() DynamicValue(
	Field<MeshType>& field,
	const BoundaryPatch& patch,
	const Time& timeArg,
	TimeSeries&& valuesArg,
	bool isUniformArg
) :
	DynamicValue( field, patch, timeArg )
{
	values( std::move(valuesArg), isUniformArg );
}

CHEAD() DynamicValue(
	Field<MeshType>& field,
	const BoundaryPatch& patch,
	const Time& timeArg,
	std::vector<FunctionParserWrapper>&& functionsArg,
	bool isUniformArg,
	const Field<MeshType>& zFieldArg
) :
	DynamicValue( field, patch, timeArg )
{
	functions( std::move(functionsArg), isUniformArg );
	zField( zFieldArg );
}

CHEAD(void) update_function() const {
	assert( isFunction() );
	assert( !functions().empty() );
	assert(
		field().nComponents() == static_cast<Index>( functions().size() )
	);
	/* a uniform function definition -> all ghost cells get the same 
		* value (which might be a vector, depending on the field) */
	if ( isUniform() ){
		for ( int i{0}; i<field().nComponents(); ++i ){
			m_uniformFunctionBuf[i] = functions()[i]
				.eval( 0, 0, 0, time().current(), time().step(), 0, 0, 0, 0 );
		}
		field().ghostValues( patch() ).colwise() = m_uniformFunctionBuf;
	/* a non-uniform function definition -> the function for each value has
	 * to be evaluated for each ghost cell */
	} else {
		/* it is valid to pass an empty z field, because not every solver
		 * uses three coordinates */
		setPatchByFunction( field(), patch(), zField(), time(), functions() );
	}
}

CHEAD(void) update_timeSeries() const {
	assert( !isFunction() );
	assert( !values().empty() );
	auto iter { findValuesAtTime( values(), time().current() ) };
	/* check whether update is needed at all. If the last value that was
	 * set is the same that would be set right now, there's no use in
	 * setting it again. */
	if (
		!time().atStart() &&
		time().current() - time().step() - iter->first > 2*epsilon
	)
		return;

	if ( isUniform() ){
		assert( iter->second.rows() == field().nComponents() &&
"Size of uniform dynamic values must be equal to number of values per cell!");
		assert( iter->second.cols() == 1 );
		
		const Index& nCells { patch().size() };
		field().ghostValues( patch() ) = iter->second.replicate(1, nCells );
	} else {
		assert( (
			iter->second.rows() == field().nComponents() &&
			iter->second.cols() == patch().size()
		) );

		field().ghostValues( patch() ) = iter->second;
	}
}

CHEAD(void) update() const {
	if ( isFunction() ){
		update_function();
	} else {
		update_timeSeries();
	}
}

CHEAD(bool) isUniform() const { return m_isUniform; }
CHEAD(void) isUniform(bool b){ m_isUniform = b; }

CHEAD(bool) isFunction() const { return m_isFunction; }
CHEAD(void) isFunction(bool b){ m_isFunction = b; }

CHEAD(auto) values() const ->
	const std::vector<std::pair<scalar,ArrayXXs>>&
{
	assert( !m_isFunction &&
		"Values unavailable when condition is set by a function!" );
	return m_values;
}

CHEAD(auto) functions() const ->
	const std::vector<FunctionParserWrapper>&
{
	assert( m_isFunction &&
		"Function parsers unavailable when condition is set by values!" );
	return m_functions;
}

CHEAD(void) values(
	TimeSeries&& values,
	bool isUniformArg
){
	assertSortedValues(values);
	isUniform( isUniformArg );
	m_isFunction = false;
	m_values = std::move(values);
}

CHEAD(void) functions(
	std::vector<FunctionParserWrapper>&& functions,
	bool isUniformArg
){
	isUniform( isUniformArg );
	if ( isUniform() )
		m_uniformFunctionBuf.resize( field().nComponents() );
	m_isFunction = true;
	m_functions = std::move(functions);
}

CHEAD(auto) time() const -> const Time& {
	assert( m_time );
	return *m_time;
}
CHEAD(void) time( const Time& timeArg ){
	m_time = make_observer<const Time>( &timeArg );
}

#undef CHEAD

#define INSTANTIATE_FOR_ALL_MESHTYPES(MeshType) \
template class DynamicValue<MeshType>;
#include "meshInstantiator.hpp"



template<typename MeshType>
void setPatchByFunction(
	Field<MeshType>& field,
	const BoundaryPatch& patch,
	const Field<MeshType>& z,
	const Time& time,
	const std::vector<FunctionParserWrapper>& functions
){
	for ( Index i{0}; i<patch.size(); ++i ){

		Index domainCell { patch.domainCells()(i) };
		auto dCentroid { field.mesh().cell(domainCell).centroid() };
		scalar zDomain { z.empty() ? 0 : z.values()(0, domainCell) };

		Index ghostCell { patch.ghostCells()(i) };
		auto gCentroid { field.mesh().cell(ghostCell).centroid() };
		/* if z is set by function, it may lead to unexpected behaviour,
		 * if the old z value is used to set the new one.
		 * In case of the first time step, the ghost z value may not be
		 * what the user intends. Therefore we check for identity.
		 * This is not a problem for zDomain, as domain values are always
		 * set first and aren't overwritten here. */
		scalar zGhost {
			(z.empty() || &z == &field) ? 0 : z.values()(0, ghostCell)
		};

		for ( Index j{0}; j<field.nComponents(); ++j ){
			field.values()( j, ghostCell ) =
				functions[j].eval(
					gCentroid.x(),
					gCentroid.y(),
					zGhost,
					time.current(),
					time.step(),
					dCentroid.x(),
					dCentroid.y(),
					zDomain,
					field.values()(j, domainCell)
				);
		}
	}
}


} // end namespace hms
